External fields, density functionals, and the Gibbs inequality 
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By combining the upper and lower bounds to the free energy as given by the Gibbs inequality for 
two systems with the same intermolecular interactions but with external fields differing from each 
other only in a finite region of space F, we show that the corresponding equilibrium densities must 
also differ from each other somewhere in F. We note that the basic equations of density functional 
theory arise naturally from a simple rearrangement and reinterpretation of the terms in the upper 
bound Gibbs inequality for such systems and briefly discuss some of the complications that occur 
when the intermolecular interactions of the two systems also differ. 
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INTRODUCTION 



Michael Fisher's work has often shown us that im- 
portant insights can arise from simple starting points 
through the use of basic principles of statistical mechan- 
ics. With this in mind, (but on a much more modest 
scale!), we will examine here some consequences of one 
of the earliest fundamental relations of statistical me- 
chanics: the Gibbs inequality, through which Gibbs es- 
tablished the basic convexity properties of the free energy 
[jlj . The usefulness of this inequality in the development 
of density functional theory (both classical ^ ||, [|, |f| 
and quantum versions 0) and in perturbation theo- 
ries of liquids J|, ||, |J has long been recognized (l(J . We 
report here some very simple consequences of the Gibbs 
inequality for the density response of a fluid to an ex- 
ternal field that do not appear to be widely known, and 
show that the basic equations of density functional theory 
arise directly from a rearrangement and reinterpretation 
of terms in a special case of the Gibbs inequality. 

The Gibbs inequality in its most general form relates 
the equilibrium properties of a system of interest with 
Hamiltonian H(<&) (with general intermolecular poten- 
tials and a generalized single particle potential $ simply 
related to the external field cj), as discussed below) to 
those of a "trial" system with Hamiltonian H t (&t) with 
possibly different intermolecular potentials and field <i> t . 
We consider here the Grand Canonical ensemble, where 
the free energy SI = — fc^TlnS and density distribution 
functions such as the singlet density p(r) are functions 
of the temperature T, volume V and chemical poten- 
tial /i, and functionals of the external field (f>(r) and of 
the (pair and any higher order) intermolecular poten- 
tials. Here 5 is the Grand partition function and we 
will let (3 — (ksT)^ 1 . We are particularly interested in 
the functional dependence on the external field 4>(r) for 
fixed values of the thermodynamic variables and the in- 
termolecular potentials. Since the zero of energy can be 
chosen arbitrarily without affecting the physics, and a 
constant field acts like a chemical potential shift in the 
Grand partition function, the relevant quantities actually 



depend only on the difference |3| 

$(r) = (f>{v) — fi 



(1) 



as shown in the Appendix. We refer to ^(r) as the sin- 
glet field and write n ([«&]) or p(r;[$]) when we want to 
emphasize this functional dependence. 

In this case the Gibbs inequality can be written as: 

Sit + (#(*) - H t (Q t )) <n<Q t + (JT($) - H t ($t))t. 

(2) 

Here ( ) and ( )t denote normalized ensemble averages 
in the systems with Hamiltonians H(&) and H t ($t) re- 
spectively, and SI and Sit are the respective free energies. 
This result follows immediately from Eq. (4) of the basic 
paper by Mermin ||; it holds for general quantum sys- 
tems and, with modified definitions of the free energies 
and averages, for other ensembles too. For complete- 
ness and to establish notation, we give in the Appendix 
a derivation for classical systems that uses only convex- 
ity properties of the exponential function and seems a 
little simpler than standard approaches ||, [|, [| [Hj. As 
emphasized in the Appendix, the inequalities in Eq. (|^) 
are strict; only if -/?($) = H t (<$>t) for every configuration 
with nonzero weight in the averages will f2 t = fi. 



II. DENSITY RESPONSE TO DIFFERENT 
SINGLET FIELDS 

Consider first the special case where the trial system 
has exactly the same intermolecular interactions as the 
system of interest and only the singlet fields may differ. 
By combining the lower and upper bounds in Eq. (Q) for 
this case we find a simple inequality 



< J dr[$(r)-$ t (r)][p t (r)-p(r)] 



(3) 



Here pt(r) = p(r; [$t]) is the equilibrium singlet density 
arising from $ t and similarly p(r) = p(r; [$]). Equation 
(^) is valid for any $t(r). Since the inequalities in Eq. (||) 
are strict, the right side of Eq. (^) is strictly greater than 
zero unless $t(r) = <&(r) for all r. 
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Equation (g) has some interesting consequences for the 
density response of a system to an external field. In par- 
ticular, it shows explicitly that the mapping from the 
singlet field $(r) to the associated equilibrium density 
p(r; [$]) is one to one: a given singlet field produces a 
unique density response. Thus if differs from <&(r) 

for any r then the right side of Eq. (ph is strictly greater 
than zero; for this to be true p*(r) must differ from p(r) 
for some r. This result plays a central role in density 
functional theory and is well known. Mermin proved this 
for quantum systems by contradiction || , and most later 
discussions of density functional theory |[ ||, || for 
classical systems have also relied on proofs by contradic- 
tionjllj. 

However Eq. (y|) allows us very simply to find more 
explicit results in some special cases that do not appear 
to be widely known. For example, suppose < &t(r) differs 
from <!>(r) only in a finite region of space T, which can 
be arbitrarily small. Then the integration in Eq. (J2J) is 
restricted to T and we find the more local result that 
Pt(r) must differ from p(r) for some r in T. (This ap- 
proach gives no information about what happens out- 
side r.) This result is obvious in the low density limit 
where there is a local density response proportional to the 
Boltzman factor of <i>(r) but it holds true in general. If in 
addition <I>t(r) > $(r) in T then J r dr[p(r) — p t (r)} > 0. 
More generally, if T denotes the set of points where $t (r) 
differs from <&(]"), then pt(r) must differ from p(r) at least 
for some r in T. 

One well known result follows when A<I> t (r) = $ t (r) — 
$(r) is very small, as indicated by the notation <5$f(r). 
Then we can use linear response theory to determine the 
small induced density change: 

Sp t (r) =-p J dr' X (v, r'; [$])5$ t (r / ), (4) 

and Eq. (j|) reduces to 

0< J drdr'5$ t (r)x(r,r';[$])5$ t (r'), (5) 

which expresses that fact that the linear response func- 
tion x is positive definite (J, [H) (except possibly at a phase 
transition in zero field where our assumption that linear 
response theory is accurate for sufficiently small S$ t (r) 
can break down). 

III. DENSITY FUNCTIONAL THEORY 

The key idea in density functional theory || ||, |], || ||, 
f?J is to consider generalized free energies that arc func- 
tionals of the density pt(v) rather than the field ^((r). 
Since the Gibbs inequality shows that a given singlet field 
produces a unique density response as discussed above, 
formally any functional of <&< can equally well be thought 
of as a functional of p t ■ But the Gibbs inequality can play 
a more central role. We will show that the basic function- 
als and the minimum principle used in density functional 



theory arise directly from the Gibbs inequality and its 
use permits us to see explicitly how the change of vari- 
able from singlet fields to densities can be carried out. 

We start from the upper bound to the free energy 
that follows from Eq. (|^) in the special case where only 
the singlet fields differ: 



il < n t + / dr[$(r)-$ t (r)]p t (r). 



(0) 



The standard density functional relations ^ |[ [|, || fol- 
low immediately by simply rearranging and reinterpret- 
ing terms in this equation. Thus we rewrite Eq. (^|) as 



Q = min 



n t - / dr Pt (r)$t(r) 



drp t (r)$(r) 



(7) 



here the minimum is taken over all possible fields <E>t and 
is achieved only when <E>t = 4> for all r. 

Equation (Q) may seem to be a particularly hard and 
inefficient way to calculate f2 but we can proceed as fol- 
lows. Since the density is related to the functional deriva- 
tive of the free energy by 



p«(r)=p(r; [* t ]) = <Jn([$ t ])/5$ t (r), 



(8) 



the terms in curly brackets in the right side of Eq. (0) 
represent a (functional) Legendre transform Q from 
n t = f2([$t]) to a so-called intrinsic free energy density 
functional F([pt\), where 



F([ Pt }) = n([$ t 



drp t (r)$ t (r) 



(9) 



denotes the terms in curly brackets. As the notation 
indicates, F is a functional of the equilibrium density 
Pt{r) and not of the field; thus the $t(r) on the right 
side is the field that corresponds to the density p*(r). 
By standard properties of the Legendre transform the 
associated field formally satisfies 



SF([p t ])/6p t (r) = -$ t (r). 



(10) 



The transformation of variables from fields to densities 
in F can be seen more explicitly if we rewrite Eq. ([)[) in 
the following equivalent form: 



F([p t }) = min 



- / dr Pt (r)<S>' t (r) 



(11) 



The minimum is taken over all possible external fields 
at fixed pt (r) . The minimum occurs when the functional 
derivative with respect to $J of the term in square brack- 
ets vanishes; from Eq. (||) this occurs when <&' t equals the 
field $ t associated with the density distribution /3t(r), as 
in Eq. (|^). Thus F(\p t ]) is certainly well-defined in prin- 
ciple, though Eq. ( |ll|) does not present a very useful way 
to calculate it in practice. 

Having made this change of variable, the only remain- 
ing dependence on the field $t on the right side of Eq. (0) 
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is through the associated equilibrium density pt(r), and 
the minimum can be taken over all such density distri- 
butions. Thus Eq. (0) can be rewritten as 



fi = 0([$] 



mm 



F([p t }) + / dr Pt (r)$(r) 



(12) 



This is the basic equation of density functional theory. 
The term in square brackets on the right side is usually 
interpreted as a generalized Grand free energy density 
functional J7$([p t ]) B: 



tlM) = F([p t }) + / drp t (r)$(r) 



(13) 



Note that ^l<s>([pt}), unlike F([p t }), is a functional of the 
field $ as well as the density p t . (Of course, both f2$ and 
F remain functionals of all the intcrmolccular interac- 
tion potentials). Equation ( |l2| ) shows explicitly that the 
minimum value of £l$([pt]) with respect to pt at fixed $ 
is the desired free energy Q ([<&]). The density that min- 
imizes this functional is p(r) and at that minimum we 
have 



0. 



(14) 



This derivation highlights the following points. The 
minimum in Eq. ( |12j ) must be taken over equilibrium 
densities pt arising from a system with the same inter- 
molecular interactions but a different singlet field. While 
the mapping from arbitrary fields to equilibrium densities 
is one to one, it is not necessarily true that an arbitrary 
function of the density must correspond to the equilib- 
rium density of a system with a given H in some field 
(unless the field is trivially taken to be so much stronger 
than the intermolecular interactions that there is an es- 
sentially local response). Thus, a density distribution 
predicting 14 nearest neighbors about a fixed particle in 
a hard sphere system is unphysical and should not be in- 
cluded in the trial density set for Eq. (|l2|). In practice, we 
suspect this is not a serious problem in most cases where 
only limited density variations are considered and for non 
hard core systems the inverse mapping can generally be 
carried out jl2| . 

Analogous equations hold in the Canonical ensemble, 
but then the trial density set should include only densities 
yielding the fixed number of particles. In this case work- 
ers have found it more convenient to work in the Grand 
ensemble where this restriction does not apply and then 
subtract out the lowest order effects of particle number 
fluctuations lO. 



IV. DIFFERENT INTERMOLECULAR 
INTERACTIONS 

More problematic is the fact that the functional de- 
pendence on p t is also affected by the intermolecular in- 
teractions. While the exact F([p t ]) is a "universal" func- 
tional of pt for a given set of intermolecular interactions, 



it remains an essentially unknown functional of the inter- 
molecular interactions, which in classical applications can 
have a variety of forms. Any approximate expression for 
F([p t ]) must implicitly or explicitly take this functional 
dependence into account. Thus a density functional that 
can accurately describes the properties of a nonuniform 
hard sphere fluid may not be appropriate for a system 
with softer repulsive intermolecular interactions. This 
problem seems much more severe when there are attrac- 
tive intermolecular interactions. Then F must describe 
critical properties, capillary waves, and a variety other 
phenomena where the functional dependence on w plays 
an essential role, and we have little idea of its form [|l4j . 
(And of course knowing the exact form, as we do for the 
original free energy 17 ([<!>]) through the partition function, 
does not mean we can carry out practical calculations!) 

In many cases where long wavelength fluctuations are 
not important, a mean field treatment of the attractive 
interactions may suffice. For the simple Lennard- Jones 
(LJ) fluid with pair interactions WLj{r 12) this can arise 
by approximating the structure of the nonuniform LJ 
fluid by that of a reference fluid with purely repulsive 
pair interactions «o(ri 2 ) giving repulsive forces equal to 
those of the LJ potential and a singlet field <E>o( r ) that 
incorporates the averaged effects of the attractive inter- 
actions p3|. The free energy can again be estimated 
using the Gibbs inequality, but now in a more general 
form than Eq. (|) or (|lf), 

since the intermolecular pair 
potentials as well as the singlet fields differ. We use the 
subscript to denote the trial system in this special case. 
Equation (|J) then becomes 

Q t -J rfrp(r)A$(r) + |J dndr 2 p {2) (r u r 2 ) Ul (r 12 ) < 
Q < Q t -J dr/9 (r)A$(r) 

+i J dridr 2 /9o 2) (ri,r 2 )wi(ri 2 ) ■ 

(15) 

Here Ui(ri 2 ) = w LJ (r 12 ) - u (ri 2 ). 

To get an upper bound to the free energy from the last 

inequality, we need an accurate expression for the pair 
(2) 

distribution function pg (ri,r 2 ) as well as po( r )- How- 
ever, since uo( r i2) is different than WLj{ri 2 ) the Gibbs 
inequality shows there is no choice of <I>o( r ) for which 
the upper bound gives the exact free energy. Neverthe- 
less, if $o( r ) is chosen so that po( r ) closely approximates 

(2) 

p(r), it seems physically plausible that p ; (ri,r 2 ) could 
also rather closely resemble // 2 )(ri,r 2 ), at least at high 
density when correlations arising from packing of the re- 
pulsive cores are most important [14 [.When this is true, 
the upper and lower bounds in Eq. (Iq) are close to each 
other and the upper bound provides a reasonable esti- 
mate for the exact free energy. The basic equations of 
perturbation theory of uniform fluids can be established 
from this perspective and accurate results for the free 
energy of the dense LJ uniform fluid are found jl5| from 
the upper bound in Eq. (|l5|). 

An upper bound for the free energy from the last in- 
equality is guaranteed only if we use accurate values for 
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(2) 

po( r ) and Pq (ri,r 2 ). In many applications of density 

(2) 

functional theory p (ri,r 2 ) is replaced by the product 
Po( r i)po( r 2) for computational simplicity S. A mini- 
mization of the approximate free energy when this ap- 
proximation is made cannot be justified by any upper 
bound principle, though it may be useful for other pur- 
poses. In a uniform fluid this approach yields a constant 
value for $o( r )j which gives a constant reference density 
equal to the LJ density, and the resulting free energy is 
only qualitatively accurate, sometimes being larger and 
sometimes smaller than the correct result H, na . 
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(A4) 
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-0Q 
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APPENDIX A 

Consider a classical system with interaction potentials 

N N 



i/(0)=^0(r l )+^^ (2) (r. i ,r, 



(Al) 



i=i 



z<3 



We explicitly denote only the dependence of H on the 
external field 4>; the intermolecular interactions are 
general and can include 3, 4, ...,N body terms. The 
Grand partition function 5 for a system with chemical 
potential p, and temperature fc^T = is given by 



S = e~ rm = Tr e-PHW+P^ = Tr e~ pH (<&) , (A2) 



where 



Tr (•) = ^[YIA 3 ^" 1 ( dr N (- 
n J 



(A3) 



or 



-pa, 



Tre 



-^ f (*t).-)3[ff(*)-fli(*t)] 



-0Q t 



(A5) 



(A6) 



If x is a random variable and (x) denotes its average 
over any normalized probability distribution we have 
e x ~( x ) > l + x—{x) for all real x. Taking averages we find 
the familiar result {e x ) > e' x ' with the strict inequality 
holding if there is any configur atio n in the average with 
x 7^ (x). Applying this to Eq. (A6) gives 



-pa 



> e 



-f3n t -0(lH(<S>)-H t (<£ t )]) t 



This yields the final result 



n < SU + (B(*)-Ht{$t)) t 



(A7) 



(A8) 



where the equality holds only when H(&) = H t (&t) for 
all configurations in the average. Swapping H(&) and 
Ht(*&t) yields the lower bound as given in Eq. (|^). By 
redefining the averages and free energies appropriately 
these results also hold in the Canonical ensemble, and 
they are valid for quantum systems as well ||. 
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